function [negLLtot,gradtot] = dreumLL2_wg(choiceFull,ppFull,tcSplit,H,K,phi,...
    PMarkov,rho,tau,theta0,tolInner,Ncpu,Nu,tcidx,chunksize)
% Calculate log-likelihood function

global gradHistory LLHistory P23maxHistory P0maxHistory


% Set parameters
alpha = theta0(end-1);
s = theta0(end);
impossible = phi == 0;
pi = [theta0(4), theta0(5), (1 - theta0(4) - theta0(5))];


% Load data from previous runs if any exists (to speed convergence)
if exist('Vbsave.mat','file') ~= 0 
    load('Vbsave.mat');    
else
    Vbsave = cell(1,Ncpu);
    for idx1 = 1:Ncpu
    	Nutp=size(tcSplit{idx1},1);
    	for idx2 = 1:tau
       		Vbsave{idx1}(:,:,:,idx2) = zeros(H+1,K,Nutp);
    	end
    end
end

% Initialize storage variables 
Pout = cell(1,Ncpu);
kappaout = cell(1,Ncpu);
Vbout = cell(1,Ncpu);
dPdvout = cell(1,Ncpu);
dvdthetaout = cell(1,Ncpu);
dzetadvout = cell(1,Ncpu);
dzetadsout = cell(1,Ncpu);
zetaout = cell(1,Ncpu);
dzetadalphaout = cell(1,Ncpu);

% Calculate log likelihood, gradient for each data chunk
parfor idx1=1:Ncpu
    
    % Load Vb to speed convergence
    tctp=tcSplit{idx1};
	temp=Vbsave{idx1};
	Vbsavetp=cell(1,tau);
	for idx2 = 1:tau 
		Vbsavetp{idx2}=temp(:,:,:,idx2);
    end
	
    % Run inner nest, save result  
	[Vb,~] = dreumInner5(H,impossible,K,tau,phi,PMarkov,rho,...
        tctp,theta0,tolInner,Vbsavetp);
    Vbout{idx1} = Vb;
    temp = zeros(H+1,K,size(tctp,1),tau);
	for idx2 = 1:tau
		temp(:,:,:,idx2) = Vb{idx2};
	end
	Vbsave{idx1} = temp;
    
    % Calculate choice probabilities, stationary distribution
    [~,kappaout{idx1},Pout{idx1},zetaout{idx1}] = dreumStationary_6(alpha,H,...
        impossible,K,phi,pi,s,tau,Vb);
    
    % Calculate gradient terms
    dPdsout{idx1} = gradPs(alpha,H,impossible,K,phi,pi,PMarkov,Pout{idx1},...
        rho,tau,tctp,theta0,tolInner,Vbsave{idx1}); % Approximate gradient wrt exponential scale parameter
    [dPdvout{idx1},dvdthetaout{idx1},dzetadalphaout{idx1},dzetadvout{idx1},dzetadsout{idx1}] ...
        = gradLL1(dPdsout{idx1},H,K,phi,rho,s,tau,tctp,theta0,Pout{idx1},...
        impossible,Vb); % Analytical gradients for all other parameters
    
end
    
% Combine terms for likelihood, gradient calculation
P = zeros(H+1,K,Nu,tau);
kappa = zeros(1,K,Nu,tau);
zeta = zeros(K,1,Nu,tau);
Vb = zeros(H+1,K,Nu,tau);
dPds = zeros(H+1,K,Nu,tau);
dzetadv = zeros(K,H-1,Nu,tau);
dPdv = zeros(H+1,H-1,K,Nu,tau);
dvdtheta = zeros(H-1,25,Nu,tau);
dzetads = zeros(K,1,Nu,tau);
dzetadalpha = zeros(K,1,Nu,tau);

Nuctr2 = 0;
for idx1 = 1:Ncpu
    n = chunksize(idx1);
    for idx3 = 1:n 
        Nuctr2 = Nuctr2+1;
        for idx2 = 1:tau
            
            P(:,:,Nuctr2,idx2) = Pout{idx1}{idx2}(:,:,idx3);
            kappa(:,:,Nuctr2,idx2) = kappaout{idx1}{idx2}(idx3,:);
            zeta(:,:,Nuctr2,idx2) = zetaout{idx1}{idx2}(:,:,idx3);
            Vb(:,:,Nuctr2,idx2) = Vbout{idx1}{idx2}(:,:,idx3);
            dPds(:,:,Nuctr2,idx2) = dPdsout{idx1}{idx2}(:,:,idx3); 
            dzetadv(:,:,Nuctr2,idx2) = dzetadvout{idx1}{idx2}(:,:,idx3);
            dPdv(:,:,:,Nuctr2,idx2) = dPdvout{idx1}{idx2}(:,:,:,idx3);
            dvdtheta(:,:,Nuctr2,idx2) = dvdthetaout{idx1}{idx2}(:,:,idx3);
            dzetads(:,:,Nuctr2,idx2) = dzetadsout{idx1}{idx2}(:,:,idx3);
            dzetadalpha(:,:,Nuctr2,idx2) = dzetadalphaout{idx1}{idx2}(:,:,idx3);
                        
        end

    end

end


% Calculate log likelihood    
[negLLtot,LL] = dreumLL2_calculate2(alpha,choiceFull,ppFull,tcidx,phi,...
    P,pi,tau,K,zeta);
    

% if nargout > 1
   
    % Calculate gradient 
    gradtot = transpose(gradLL2_2(alpha,choiceFull,dPds,dPdv,dvdtheta,dzetadalpha,...
        dzetads,dzetadv,H,K,LL,phi,ppFull,tau,tcidx,P,pi,zeta));   
    
% end

gradHistory = [gradHistory; gradtot'];
LLHistory = [LLHistory; -negLLtot];
P23maxHistory = [P23maxHistory; max(max(P(end,:,:,:)))];
P0maxHistory = [P0maxHistory; max(max(P(1,:,:,:)))];
save('gradHistory.mat','gradHistory');
save('LLHistory.mat','LLHistory');
save('P23maxHistory.mat','P23maxHistory');
save('P0maxHistory.mat','P0maxHistory');

% % Approximate gradient (delete when running for real)
% (negLLtot-negLLtotKEEP)/diff

save('Vbsave.mat','Vbsave')